Statistical distributions¶
Fermi-Dirac statistics, Bose-Einstein statistics, Maxwell-Boltzmann distribution, Planck’s law Binomial distribution
- number of successes of n independent experiments Cauchy distribution
- ration of normally distributed variables
- shape of spectral lines Gamma distribution Geometric distribution
- number of trials before the first success in Bernoulli trials Chi-squared distribution
- sum of squares of independent standard normal variables
- how to construct confidence intervals
Exponential distribution Gumbel distribution
- expectation value of the maximum/minimum of samples drawn from a parent distribution with exp tail
- Lévy distribution
- heavy-tailed distribution
Logistic distribution
- categorical dependent variables
- logistic regression
Log normal distribution
- product of many independent variables
- size distribution of particles
Normal distribution
Pareto distribution
- heavy-tailed, decreasing power-law
- degree distribution of social networks
Poisson distribution
Student t
- average of a normally distributed population
- n is small
Uniform distribution
Weibull distribution
- any of n iid variables being above a certain limit
- first failure of a system of n identical components
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import plotly.graph_objects as go
import pandas as pd
import seaborn as sns
from copy import deepcopy
from bs4 import BeautifulSoup
from urllib.request import Request, urlopen
from scipy.stats import gumbel_r, gompertz, maxwell, pareto, levy, lognorm, weibull_max, expon, poisson
matplotlib.rcParams["figure.figsize"] = (10, 6)
matplotlib.rcParams["axes.titlesize"] = 16
matplotlib.rcParams["axes.labelsize"] = 13
matplotlib.rcParams["xtick.labelsize"] = 16
matplotlib.rcParams["ytick.labelsize"] = 16
matplotlib.rcParams["axes.grid"] = True
matplotlib.rcParams["xtick.major.size"] = 12
matplotlib.rcParams["ytick.major.size"] = 12
matplotlib.rcParams["xtick.major.width"] = 2
matplotlib.rcParams["ytick.major.width"] = 2
Excercise 1¶
For the purposes of the current excercises, select a stream from a rainy area with relatively small discharge so that the effect of short but strong storms is visible. Good choices are small rivers from the north-eastern US, e.g. site 01589440. Retrieve at least 10 years of data.
Large amounts of historical surface water data are available from the United States Geological Survey (USGS) at https://waterdata.usgs.gov/nwis The goal of the project is to retrieve samples from the web interface manually, and then later automate the process by calling the web service as described at https://help.waterdata.usgs.gov/faq/automated-retrievals.
You can also use the dataset from the /home/course/Datasets/02-data/ folder as well.
or¶
Find a database of historic water level measurement of the Danube at Budapest
my_url = "https://nwis.waterdata.usgs.gov/usa/nwis/uv/?cb_00060=on&cb_00065=on&format=html&site_no=03014500&period=&begin_date=2010-02-08&end_date=2021-02-15"
#Downloading the webpage and storing its HTML right away.
resp_obj = urlopen(Request(my_url, headers={'User-Agent': 'Mozilla/6.0'}))
page_html = resp_obj.read()
resp_obj.close()
#Parsing the website's HTML with bs4. For some reason, I couldn't reach the body.table part in the traditional way.
#So, I had to go with a bruteforce method. Of course, if the table would've been consisted of more columns with
#other types of data, this code wouldn't have worked right away.
parsed_text = BeautifulSoup(page_html,"html.parser")
rows = [x for x in parsed_text.body.contents[11].contents[4].children]
rows = list(filter(lambda x: x != "\n", rows)) #filtering out newline characters
data = {} #Storing every row of the table in a key-value pair (dict)
for i in range(len(rows)):
row = rows[i]
datetime_timezone = row.findAll("td", {"nowrap":"nowrap"})[0].string
discharge_value = row.findAll("td", {"nowrap":"nowrap"})[1].span.string
discharge_comment = row.findAll("td", {"nowrap":"nowrap"})[1].sup.string
gage_height_value = row.findAll("td", {"nowrap":"nowrap"})[2].span.string
gage_height_comment = row.findAll("td", {"nowrap":"nowrap"})[2].sup.string
data[f"row_{i}"] = [datetime_timezone, discharge_value, discharge_comment, gage_height_value, gage_height_comment]
data_df = pd.DataFrame.from_dict(data, orient="index", columns = ["datetime_tz", "discharge_v", "discharge_c", "gage_h", "gage_c"])
data_df.head()
Automated retrieval¶
!curl "https://nwis.waterdata.usgs.gov/nwis/uv/?parameterCd=00060,00065&format=rdb&site_no=01321000&period=&begin_date=20100101&end_date=20191231&siteStatus=all" > data.csv
!head -n 40 data.csv
Excercise 2¶
Load the downloaded data file into the processing environment paying attention to handling time stamps and perfoming the necessary data type conversions. Converting dates to floating point numbers such as unix time stamp or julian date usually makes handling time series easier. Plot the data for a certain interval to show that the effect of storms is clearly visible.
df = pd.read_csv("data.csv",sep='\t', comment='#',header=0)
df = df.drop(index=0)
df.head()
print(df.agency_cd.unique(), df.site_no.unique(), df.tz_cd.unique())
#These columns are irrelevant, because they contain the same data
df = df.drop(columns = ["agency_cd", "site_no","tz_cd"])
df.head()
df.dtypes #The columns has to be cast into the appropriate data formats
From the head of the .csv datafile we know what the codes are
106704 00065 Gage height, feet
106705 00060 Discharge, cubic feet per second
df = df.rename(columns = {"106704_00065": "gage", "106704_00065_cd": "gage_c", "106705_00060": "discharge", "106705_00060_cd":"discharge_c"})
df["datetime"] = pd.to_datetime(df.datetime)
df["discharge"] = df.discharge.astype(float)
df["gage"] = df.gage.astype(float)
df.head(), df.dtypes
#We can see that there are some NaN values, the ratio of these NaNs in each columns are the following:
df.isna().sum()/len(df)
df = df[df.discharge.isna() != True] #Getting rid of the rows where discharge is missing
df.isna().sum()/len(df)
Excercise 3¶
Plot the histogram of water discharge values. Fit the data with an appropriate distribution function and bring arguments in favor of the choice of function.
def fitPlotDist(dist,to_fit,x, alpha = 1.0):
params = dist.fit(to_fit)
rv = dist(*params)
plt.plot(x, rv.pdf(x), label=dist.name, alpha = alpha)
sns.distplot(df.discharge,bins=2000, label="raw data with KDE", kde = True)
plt.ylabel("frequency")
x = np.linspace(min(df.discharge), max(df.discharge), 1000)
#fitPlotDist(gumbel_r,df.discharge,x)
#fitPlotDist(pareto, df.discharge, x)
#fitPlotDist(maxwell, df.discharge,x)
fitPlotDist(lognorm,df.discharge,x)
fitPlotDist(levy, df.discharge, x)
fitPlotDist(expon,df.discharge,x)
plt.xlim(0,5000)
plt.legend()
plt.title("Histogram of water discharge values")
The data visualized above shows, that from the three distribution I've fitted it, the Lévy distribution is the best fit. I didn't find anything about the Lévy distribution on Wikipedia that would support this claim.
However, I suspect, that the lognormal distribution would be more well-founded, because as Wikipedia states:
a log-normal process is the statistical realization of the multiplicative product of many independent random variables.
Which is true in the case of discharge values: it is affected by many variables be it meteorological or geological.
Excercise 4¶
In case of small streams, storms and passing weather fronts with rain can significantly increase water discharge on a short time scale. Develop a simple algorithm to detect rainy events in the data. Plot the original time series and mark rainy events.
fig, ax1 = plt.subplots(figsize=[15,5])
ax1.plot(df.datetime,df.discharge, color="tab:orange", label="Discharge")
ax1.set_ylabel("Discharge [$ft^3$/s]",size=13)
plt.grid(ls = "dashed", color="tab:orange")
plt.legend(loc = 2,fontsize=12)
ax2 = ax1.twinx()
ax2.plot(df.datetime, df.gage, color="tab:blue",alpha = 0.6, label="Gage height")
ax2.set_ylabel("Gage height [ft]", size=13)
plt.grid(ls = "dotted",alpha = 0.5,color="tab:blue")
plt.legend(loc = 1, fontsize=12)
plt.title("Discharge and gage height plotted against time")
from scipy.stats import pearsonr
pearsonr(df[df.isna().sum(axis=1) < 1].discharge, df[df.isna().sum(axis=1) < 1].gage)
From this graph we can conclude, that the discharge at a given cross section is linearly dependant with the gage height, which seems intuitive. This claim is also supported by the strong positive correlation score.
'''fig, ax1= plt.subplots(figsize=[15,5])
ax1.plot(df.datetime,df.discharge.diff(), color="tab:orange", label="Change in discharge")
ax1.set_ylabel("Change in discharge [$ft^3$/s]",size=13)
plt.grid(ls = "dashed", color="tab:orange")
plt.legend(loc = 2)
ax2 = ax1.twinx()
ax2.plot(df.datetime, df.gage.diff(), color="tab:blue",alpha = 0.6, label="Change in gage height")
ax2.set_ylabel("Change in gage height [ft]", size=13)
plt.grid(ls = "dotted",alpha = 0.5,color="tab:blue")
plt.legend(loc = 1)
plt.title("Change in discharge and in gage height plotted against time",size=16)'''
from scipy.signal import argrelextrema, find_peaks
plt.figure(figsize=[15,5])
plt.plot(df.datetime,df.discharge, color="tab:orange", label="Discharge")
#p = argrelextrema(df.discharge.to_numpy(), np.greater, order = 5)[0]
p,_ = find_peaks(df.discharge,threshold=None, height = 2000, prominence = 700)
plt.plot(df.datetime.to_numpy()[p],df.discharge.to_numpy()[p], "o", label="Possible rain",color="tab:green")
plt.ylabel("Discharge [$ft^3$ /s ]",size=13)
plt.legend()
#plt.ylim(0,3000)
plt.title("Possible rains based on sudden increase in discharge",size=16)
len(p)
fig = go.Figure()
fig.add_trace(go.Scatter(x=df.datetime,
y=df.discharge,
name = "Discharge",
line={ "color" : "#FF7F0E", "width": 2}))
p,_ = find_peaks(df.discharge,threshold=None, height = 500, prominence = 500)
#Plot only the points, where peaks were found
fig.add_trace(go.Scatter(x = pd.to_datetime(df.datetime.to_numpy()[p]),
y = df.discharge.to_numpy()[p],
mode='markers',
name='Possible rain',
marker = { "color" : "#2CA02C"}))
fig.update_layout(yaxis_title="Discharge [ft^3 /s ]")
#fig.update_xaxes(range=[pd.to_datetime("2014-01-01 00:00:00"), pd.to_datetime("2015-12-31 23:45:00")])
#This was my test interval
fig.show()